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PREFACE 


This  report  is  sponsored  by  the  Office,  Chief  of  Engineers,  U.  S. 
Amy,  OCE,  as  part  of  the  Environaental  Water  Quality  and  Operational 
Studies  (EWQOS)  Work  Unit  IC.l,  entitled  Improve  and  Verify  Existing 
One-Dimensional  Reservoir  Water  Quality  and  Ecological  Predictive 
Techniques.  The  OCE  Technical  Monitors  for  EWQOS  were  Mr.  Earl  Eiker, 
Mr.  John  Bushman,  and  Mr.  James  L.  Gottesman. 

Work  for  this  report  was  conducted  during  the  period  October  1981- 
March  1982  by  Dr.  Joseph  H.  Wlosinski,  Water  Quality  Modeling  Group 
(WQMG)  of  the  Environmental  Laboratory  (EL).  The  draft  report  was  re¬ 
viewed  by  Drs.  Carol  Collins  and  Allan  S.  Lessem,  Mr.  Jack  B.  Waide, 
and  Ms.  Linda  S.  Johnson. 

The  study  was  conducted  under  the  direct  supervision  of  Mr.  Aaron 
Stein,  Acting  Chief,  WQMG,  and  Mr.  D.  L.  Robey,  Chief,  Ecosystem  Re¬ 
search  and  Simulation  Division,  and  under  the  general  supervision  of 
Dr.  John  Harrison,  Chief,  EL,  U.  S.  Army  Engineer  Waterways  Experiment 
Station  (WES).  Program  Manager  of  EWQOS  was  Dr.  J.  L.  Mahloch,  EL. 

Director  of  WES  during  this  study  and  the  preparation  of  this 
report  was  COL  Tilford  C.  Creel,  CE.  Technical  Director  was 
Mr.  F.  R.  Brown. 


This  report  should  be  cited  as  follows: 


Wlosinski,  J.  H.  1984.  "Evaluation  Techniques  for 
CE-QUAL-R1:  A  One-Dimensional  Reservoir  Water  Quality 
Model,"  Miscellaneous  Paper  E-84-1,  u.  S.  Army  Engineer 
Waterways  Experiment  Station,  Vicksburg,  Miss. 
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EVALUATION  TECHNIQUES  FOR  CE-QUAL-R1 :  A  ONE-DIHENSIONAL 
RESERVOIR  WATER  QUALITY  MODEL 

PART  I :  INTRODUCTION 

Background 

1.  One  of  the  highest  priority  needs  of  U.  S.  Army  Corps  of 
Engineers  (CE)  District  and  Division  Offices  is  the  ability  to  realis¬ 
tically  predict  and  assess  the  effects  of  engineering  activities  on 
the  environment  (Keeley  et  al.  1978).  To  help  assess  engineering 
effects  on  reservoir  water  quality,  a  one-dimensional  mathematical  model 
called  CE-QUAL-R1  (Environmental  Laboratory  1982)  which  includes  phys¬ 
ical,  chemical,  and  biological  factors  is  being  developed  as  part  of  the 
Environmental  Water  Quality  and  Operational  Studies  (EWQOS)  Program. 

The  forerunner  of  CE-QUAL-RI  was  the  reservoir  portion  of  a  model 
called  Water  Quality  for  River-Reservoir  Systems  (WQRRS),  which  was 
assembled  in  1974  for  the  U.  S.  Army  Engineer  Hydrologic  Engineer¬ 
ing  Center,  Davis,  California,  by  Water  Resources  Engineers,  Inc. 
CE-QUAL-RI,  which  resides  on  the  Boeing  Hainstream-EKS  interactive 
time-sharing  computer  system,  has  been  used  in  the  past  to  evaluate 

p re impoundment  water  quality  problems  and  the  effects  of  reservoir 
operation  on  water  quality  (see,  for  example,  Ford  et  al.  1977,  1979 
and  Thornton  et  al.  1976,  1977). 

2.  Another  task  within  the  EWQOS  Program  includes  long-term 
comprehensive  reservoir  field  studies  (Work  Unit  VIIA).  These  field 
studies  have  provided  data  which  are  especially  suitable  for  use  with 
CE-QUAL-RI.  Data  from  DeGray  Reservoir,  Arkansas,  and  Eau  Galle  Reser¬ 
voir,  Wisconsin,  have  been  used  here  to  provide  information  for  improv¬ 
ing  predictive  capabilities  for  CE-QUAL-RI;  these  data  were  collected 
at  biweekly  or  monthly  intervals,  usually  at  meter  increments  of  depth, 
a  scheme  which  is  suitable  for  evaluating  the  model. 
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Purpose 


3.  The  purpose  of  this  report  is  to  discuss  methods  to  be  used 
for  evaluating  the  mathematical  model  CE-QUAL-R1.  Tests  are  proposed 
to  ensure  that  the  coding  of  the  model  is  correct  and  to  ensure  that 
model  predictions  are  suitable  for  the  needs  of  CE  District  and  Divi¬ 
sion  Offices.  The  proposal  includes  comparisons  of  model  predictions 
with  field-measured  values.  The  actual  utilization  of  these  methods  to 
evaluate  the  application  of  CE-QUAL-R1  to  DeGray  and  Eau  Galle  Reser¬ 
voirs  will  be  discussed  in  subsequent  reports. 
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PART  II:  LITERATURE  REVIEW 


4.  Although  the  literature  contains  much  information  concerning 
model  evaluation,  there  is  little  agreement  among  authors.  In  fact, 
most  modelers  cannot  agree  on  terminology.  Not  only  are  different  terms 
used  for  the  same  process,  but  the  same  term  is  used  for  different  pro¬ 
cesses.  For  example,  the  process  of  comparing  model  output  to  field- 
measured  values  is  referred  to  as  validation  by  a  number  of  authors 
(e.g. ,  House  1974,  Miller  1974,  Hall  and  Day  1977,  Schruben  1980,  and 
Gentil  and  Blake  1981)  but  as  verification  by  others  (e.g.,  Orlob  1975, 
Weatherbe  1976,  Bedford  and  Baba j imopoulos  1980,  Thomann  1980,  and 
Reckhow  1981).  In  addition,  verification  has  been  used  by  House  (1974) 
and  Mihram  (1972)  to  describe  a  means  to  test  the  consistency  of  model 
design  or  its  intended  algorithmic  structure,  whereas  the  same  general 
process  is  termed  validation  by  Lawler  (1980).  Goodall  (1972)  also 
used  the  term  validation,  but  he  suggested  that  the  process  will  not 
tell  us  if  a  model  is  valid  or  invalid.  O'Neill  (1975)  wrote  that  it 

is  possible  to  invalidate  or  validate  the  same  model  by  manipulation 
of  the  questions  asked.  Caswell  (1976)  wrote  that  predictive  models 
should  be  validated,  and  theoretical  models  corroborated .  In  reviewing 
Caswell's  paper,  Wiegert  (1975)  accepted  the  term  corroboration  but 
said  that  acceptance  might  be  a  more  preferable  term  than  validation. 
Nolan  (1972)  verified  coding  and  assumptions  but  validated  hypotheses 
and  recognition  of  perception  filters.  Mankin  et  al.  (1977)  suggested 
that  one  should  dismiss  the  question  of  model  validity  and  ask  instead 
whether  or  not  the  model  is  useful. 

5.  Although  the  terminology,  number  of  steps,  and  methods  may 
not  be  agreed  upon  by  authors,  tests  for  the  evaluation  of  models  gen¬ 
erally  involve  two  processes.  The  first  process  tests  whether  or  not 
the  model  responds  in  the  manner  that  the  modeler  intended;  this  process 
involves  "debugging"  the  model,  but  should  include  other  tests  as  well 
(Mihram  1972).  Mihram  suggested  a  systematic  test  that  determines  some 
specific  set  of  environmental  conditions  for  which  the  model's  response 
should  be  known.  A  similar  test  was  suggested  by  Lawler  (1980),  who 
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recommended  using  repetitive  input  and  boundary  conditions  and  checking 
the  values  of  the  state  variables  after  sufficient  time  had  elapsed  for 
the  model  to  dampen  transient  behavior. 

6.  The  second  process  compares  model  predictions  with  field- 
measured  values  in  an  attempt  to  demonstrate  that  the  model  acceptably 
simulates  the  real  world.  It  is  this  process  about  which  most  has  been 
written,  but  upon  which  most  disagreement  remains.  Even  though  the 
questions  "Has  the  model  been  verified?"  or  "Has  the  model  been  vali¬ 
dated?"  are  often  asked,  a  number  of  authors  state  that  a  valid  model 
is  impossible.  Goodall  (1972)  stated  that  since  models  are  simplifica¬ 
tions,  they  will  virtually  never  be  exact  representations  of  conditions 
in  the  real  world.  He  argued  that  the  question  of  whether  or  not  a 
model  should  be  accepted  or  rejected  is  not  appropriate;  one  should  seek 
rather  to  evaluate  the  "goodness"  of  the  predictions  and  the  errors 
associated  with  the  model.  Schruben  (1980)  also  stated  that  the  devel¬ 
opment  of  a  strictly  valid  model  of  a  nontheoretical  process  is  impossi¬ 
ble.  Reckhow  (1981),  in  a  philosphical  analysis  of  model  verification, 
stated  that  the  testing  of  models  is  an  inductive  process  and  that  veri¬ 
fication,  which  he  defined  as  the  ascertainment  of  truth,  is  inconsis¬ 
tent  with  the  inductive  logic  of  scientific  research.  Wright  (1972), 

in  a  review  of  ecosystem  models,  agreed  that  predictive  models  could 
not  be  invalidated.  House  (1974)  hypothesized  that  the  world  is  so 
complex  that  attempts  to  completely  validate  forecasting  models  are 
futile. 

7.  Regardless  of  whether  or  not  one  can  state  that  a  model  is, 
or  is  not,  valid,  methods  have  been  put  forth  to  compare  model  output 
to  measured  values.  Graphical  techniques  are  often  used;  i.e.,  a  par¬ 
ticular  variable  is  plotted  through  time  as  predicted  and  as  measured, 
and  the  plots  are  compared  qualitatively.  In  addition  to  qualitative 
tests,  a  number  of  quantitative  tests  have  been  proposed.  Some  of  the 
tests  appear  to  have  been  developed  especially  for  model  testing;  these 
include  Theil's  (1961)  inequality  coefficient,  Kapoor's  (1968)  inequal¬ 
ity  coefficient,  a  reliability  index  for  models  (Leggett  and  Williams 
1981),  'nd  a  procedure  for  simulation  model  acceptance  (Schruben  1980). 
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Other  quantitative  tests  of  model  validity  are  based  on  statistical 
measures  which  generally  show  whether  or  not  two  samples  came  from  the 
same  population.  (Table  1  lists  some  of  these.)  Some  of  the  statis¬ 
tical  tests  are  appropriate  for  distributions  of  residuals,  others  for 
comparison  of  frequency  distributions;  some  are  for  deterministic 
models,  and  others  for  stochastic  models;  some  of  the  tests  are  para¬ 
metric,  while  others  are  nonparametric.  As  with  terminology,  there 
does  not  seem  to  be  general  agreement  on  particular  tests,  and  tests 
recommended  by  some  authors  may  be  rejected  by  others.  For  example, 
according  to  Wright  (1972),  "Using  either  regression  or  factor  analysis 
to  validate  computer  simulation  models  is  absurd";  concerning  tests  of 
selected  measures  from  distributions,  he  wrote  "At  worst,  such  tests  are 
totally  improper;  at  best,  they  destroy  the  structure-in-time  of  the 
trajectories  being  contrasted." 

Table  ] 

A  Partial  List  of  Statistical  Tests  for  Comparing  Model 
Predictions  and  Measured  Data* 


Analysis  of  variance 
Chi-square 
Comparison  of  means 
Factor  analysis 
F-test 
Kendall  tau 

Kolmogorov-Smirnov  test 

Mann-Whitney  test 

Normalized  mean  error 

Pearson  product  moment  cor¬ 
relation  coefficient 


Regression  analysis 
Relative  error 
Root-roean-square  error 

1  sample  t-test 

2  sample  t-test 
Sign  test 
Spearman  rho 
Spectra]  analysis 
Wald-Wolfowitz  test 
Wilcoxon  test 


*  From  Mihram  (1972),  Reckhow  (1981),  Wright  (1972),  Thomann  (1980), 
Gordon  (1981),  and  TRC  (1981). 
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8.  A  workshop  on  verification  of  water  quality  models  was  con¬ 
vened  by  the  Environmental  Protection  Agency  in  1979.  The  proceedings 
from  the  workshop  (U.  S.  Environmental  Protection  Agency  (EPA)  1980)  are 
especially  apropos  to  the  present  model  evaluation.  The  members  of  the 
workshop  recommended  the  review  of  software  code  and  the  use  of  internal 
automated  checks  for  evaluating  computer  programs.  They  also  generally 
agreed  that  an  adequate  model  evaluation  would  consist  of  comparing  com¬ 
puted  model  results  to  a  set  of  water  quality  data  other  than  the  cali¬ 
bration  data  set.  Although  they  encouraged  the  use  of  statistical  tech¬ 
niques,  they  did  not  recommend  any  statistics,  nor  did  they  believe 


PART  III:  RATIONALE  AND  EVALUATION  METHODS 


Model  Objectives 

9.  Innis  (1975)  and  Swartzman  (1979)  argued  that  evaluation  cri¬ 
teria  should  depend  on  model  objectives.  In  agreement  with  their  argu¬ 
ment,  the  objective  of  developing  the  CE-QUAL-R1  model  may  be  stated 
thus:  to  provide  CE  District  and  Division  offices  with  a  means  of  study¬ 
ing  preimpoundment  and  postimpoundment  water  quality  problems  and  the 
effects  of  reservoir  operation  on  water  quality.  Examples  of  functions 
that  the  model  can  perform  include: 

a.  Determine  onset,  extent,  and  duration  of  thermal 
stratification. 

b.  Locate  selective  withdrawal  ports  required  to  meet  a 
downstream  water  quality  objective. 

c.  Determine  the  effect  of  structural  modifications  on  water 
quality. 

d.  Predict  the  development  of  anoxic  conditions. 

e.  Provide  information  concerning  algal  blooms. 

f.  Isolate  factors  limiting  algal  growth. 

g.  Predict  effects  of  storm  events  on  inpool  and  release 
water  quality. 

h.  Determine  effects  of  upstream  land  use  on  inpool  and 
release  water  quality. 

i^.  Determine  effects  of  project  operation  changes  such  as: 

(1)  Altered  release  level. 

(2)  Change  in  minimum  or  maximum  release  rate. 

(3)  Changes  in  pool  elevation. 

(4)  Destratification. 

10.  CE-QUAL-R1  will  be  used  as  a  management  tool,  often  on  reser¬ 
voirs  that  are  only  in  the  planning  stage.  The  model  must  therefore  be 
general  enough  to  allow  for  its  use  in  simulating  a  variety  of  impound¬ 
ments,  planned  or  in  existence,  with  a  host  of  possible  operational 
plans.  CE-QUAL-R1  fulfills  that  requirement  because  it  is  not  a  model 
per  se,  using  the  same  set  of  initial  conditions,  coefficients,  and 
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updates  for  all  possible  reservoirs.  It  is  rather  a  model  framework 
which  becomes  a  model  for  a  particular  waterbody  after  it  has  incorpo¬ 
rated  the  initial  conditions,  coefficients,  and  descriptions  that  are 
characteristic  of  a  particular  site. 

Evaluation  Objective 

11.  With  the  above  two  model  objectives  in  mind,  the  goal  of  the 
evaluation  process  will  be  to  supply  the  best  possible  tool,  within  the 
assumptions  specified  for  CE-QUAL-R1,  for  reservoir  water  quality  man¬ 
agement.  The  process  will  not,  as  has  been  argued  by  Goodall  (1972), 
Schruben  (1980),  Reckhow  (1981),  Wiegert  (1975),  House  (1974),  and 
others,  tell  if  CE-QUAL-R1  is  or  is  not  valid.  There  are  unresolved 
questions  about  the  validity  of  the  predictions  of  any  model  that  pre¬ 
dicts  numerous  variables  in  a  number  of  layers  and  whose  predictions 
may  not  fall  within  some  confidence  band.  Suppose,  for  example,  pre¬ 
dictions  of  oxygen  for  1  year  were  compared  to  monthly  observed  data 
and  all  comparisons  were  satisfactory  except  those  predicted  for  May. 

If  one  considers  the  predictions  to  be  satisfactory  11  out  of  12  times, 
one  assumes  that  there  is  some  sort  of  self-correcting  mechanism  built 
into  the  model  concerning  oxygen  prediction,  for  the  predicted  value  at 
one  point  in  time  depends  on  the  previously  predicted  values.  Would 
the  predictions  be  considered  valid  only  through  April? 

12.  Suppose,  in  another  case,  that  a  model  predicts  both  oxygen 
and  algae.  Suppose  further,  that  comparisons  show  that  oxygen  predic¬ 
tions  are  satisfactory  but  algal  predictions  are  not.  Can  the  model  be 
considered  okay  for  predicting  oxygen?  Because  the  oxygen  predictions 
are  in  part  based  on  concentrations  of  algae  which  are  not  satisfactory, 
the  model  may  be  predicting  the  correct  oxygen  values  for  the  wrong 
reason. 

13.  With  the  above  arguments  in  mind,  and  in  agreement  with  the 
authors  cited  above,  this  author  does  not  believe  that  the  evaluation 
process  can  supply  a  model  that  has  been  completely  verified,  or  one 
that  is  or  is  not  valid.  He  offers,  instead,  an  evaluation  of  CE- 
QUAL-R1  that  will  improve  predictions  by  (a)  testing  alternate 
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algorithms  for  the  same  process,  or  (b)  testing  the  results  achieved 
using  alternate  processes  or  variables.  The  results  of  the  evaluation 
must  then  be  transferred  to  the  managers  who  intend  to  use  model  predic¬ 
tions,  because  the  ultimate  judgement  of  a  model  depends  on  the  objec¬ 
tives  of  a  particular  study.  As  with  most  other  model  evaluations  (see 
paragraphs  5  and  6),  two  main  processes  will  be  used  for  the  evaluation: 
the  first  tests  the  software  code;  the  second  examines  model  predictions. 

14.  The  term  calibration  will  be  used  in  this  report  to  mean  the 
process  of  adjusting  a  set  of  coefficients  by  comparing  model  predic¬ 
tions  to  measured  values  for  a  data  set  representing  a  particular  reser¬ 
voir  for  a  particular  period  of  time.  The  term  verification  will  be 
used  to  mean  the  process  of  comparing  model  predictions  to  measured  val¬ 
ues  using  a  data  set  representing  the  same  reservoir  used  for  calibra¬ 
tion,  but  for  a  different  time  period.  The  verification  data  set  must 
retain  the  same  coefficients  as  were  used  for  the  calibration  exercise. 

Evaluation  of  Software 

15.  The  origin  of  CE-QUAL-R1  dates  back  to  1972  (Environmental 
Laboratory  1982).  Since  that  time  numerous  changes  have  been  made  to 
the  code.  Due  to  the  size  of  the  model  and  the  number  of  interactions 
among  variables,  some  of  the  changes  may  have  inadvertently  caused  prob¬ 
lems  in  other  parts  of  the  model.  Some  errors  may  be  very  difficult  to 
find  because  all  of  the  code  is  not  necessarily  executed  during  every 
simulation.  Among  oOier.s,  the  following  methods  are  proposed  for  soft¬ 
ware  evaluation. 

a.  Check  the  equations  for  correct  dimensional ity . 

b.  Ensure  that  model  predictions  are  numerically  stable. 
Predictions  should  not  oscillate  more  than  measureaients 
found  in  nature  from  one  tisie  step  to  the  next;  in  addi¬ 
tion,  predictions  should  not  vary  appreciably  when  the 
tim^  step  is  varied. 

c.  Test  for  conservation  of  suss.  The  test  should  be  made 
for  both  conservative  and  nonconservative  substances. 

d.  Check  initial  values  for  zero  entries.  Occasionally  a 
variable  may  not  be  given  a  correct  initial  value,  in 
which  case  the  computer  supplies  a  zero  value.  The  zero 
value  may  allow  computations  concerning  the  variable  to 
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be  carried  out,  but  the  results  are  incorrect.  A  method 
supplied  for  the  Boeing  Mainstream-EKS  interactive  time¬ 
sharing  computer  system  allows  the  checking  of  all  vari¬ 
ables  to  ensure  that  initial  values  were  set. 

e.  As  was  suggested  by  Mihram  (1972)  and  Lawler  (1980), 
constant  values  for  a  conservative  substance  can  be  used 
for  inflowing  concentrations,  which  should  force  pre¬ 
dicted  values  in  the  water  column  to  approach  these  con¬ 
stant  values.  Increased  flows  should  cause  the  constant 
values  to  be  reached  more  quickly.  Care  must  be  taken  as 
this  test  can  apply  to  the  entire  water  column  only  dur¬ 
ing  isothermal  conditions.  The  mixing  coefficients  can 
be  changed  to  ensure  complete  mixing. 

f.  Thoroughly  check  problems  reported  by  individuals  or 
groups  which  are  using  CE-QUAL-R1 . 


Evaluation  of  Model  Predictions 

Graphical  comparisons 

16.  Models  are  here  evaluated  in  both  a  qualitative  and  quanti¬ 
tative  fashion.  Qualitatively,  the  interactive  graphics  package 
(Environmental  Laboratory  1982)  is  used,  with  predicted  and  measured 
values  graphed  together.  Figure  1  is  an  example  of  graphical  output. 


Figure  1.  An  example  of  graphical  output  comparing 
swasured  and  predicted  values  (DeCray  Reservoir 
oxygen  profile  on  1  Nay  1979) 
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Graphical  comparisons  are  used  in  addition  to  statistical  analysis 
because  statistics  can  often  be  misleading.  Consider,  for  example, 
Figure  2.  In  both  figures  the  solid  circles  represent  measured  values 
and  the  line  represents  model  predictions.  The  lines  in  the  two  figures 
represent  different  algorithms  for  the  same  process.  Even  though  the 
comparisons  in  Figure  2a  appear  to  be  better,  most  statistical  analyses 
would  show  that  the  algorithm  represented  in  Figure  2b  is  superior.  The 
reason  for  this  apparent  anomaly  is  that,  although  the  measured  points 
in  Figure  2a  are  close  to  the  predicted  line  horizontally  along  the  time 
axis,  they  are  not  so  close  as  in  Figure  2b  when  measured  vertically 
along  the  concentration  axis  at  common  points  in  time.  Since  most 
statistical  analyses  compare  the  concentrations  of  measured  and  pre¬ 
dicted  values  at  the  same  points  in  time,  the  line  in  Figure  2b  actually 
is  statistically  a  better  representation  of  measured  values. 

Statistical  analyses 

17.  Although  conclusions  about  the  adequacy  of  model  predictions 
for  a  particular  variable  can  be  made  when  viewing  a  graph,  the  number 
of  variables  and  the  number  of  layers  for  which  measured  data  are  avail¬ 
able  could  be  so  great  as  to  preclude  adequate  judgement  of  the  total 
model.  For  these  cases  statistical  analyses  would  be  useful.  Statisti¬ 
cal  packages  can  be  used  to  test  which  of  two  algorithms  for  a  particu¬ 
lar  process  is  a  better  predictor  or  which  of  a  number  of  sets  of  coef¬ 
ficients  produce  simulation  curves  which  most  closely  correspond  to 
observed  data.  A  statistics  program  is  available  for  evaluating 
CE-QUAL-R1.  It  contains  the  following  statistics:  (a)  reliability  in¬ 
dex  (Leggett  and  Williams  1981);  (b)  paired  t-test  for  means  (Sokal  and 
Rohlf  1969);  (c)  normalized  mean  error  (Gordon  1981,  see  also  Wlosinski 
1982);  and  (d)  coefficients  for  the  linear  regression  equation  for  plot¬ 
ting  observed  versus  predicted  values  (Thomann  1980).  Equations  for 
these  statistics  can  be  found  in  Appendix  A.  Note  that  these  statisti¬ 
cal  equations  can  be  applied  either  to  data  collected  at  a  single  point 
in  time,  or  to  data  collected  at  a  number  of  times  throughout  a  sam¬ 
pling  period.  In  this  way  the  evaluation  of  predicted  time  series  or 


trends  can  be  accomplished. 


CONCENTRATION 


TIME 


a.  Algorithm  1  results  compared  to  measured  values 


b.  Algorithm  2  results  compared  to  measured  values 

Figure  2.  Example  of  a  comparison  of  results  of  two  models 
to  the  same  data  set.  The  solid  line  represents  predicted 
values;  the  circles  represent  measured  values 
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IB.  Results  from  statistical  analysis  can  be  misleading,  so  the 
user  is  cautioned  to  view  graphs  of  simulation  results  in  addition  to 
using  the  statistical  package.  To  help  familiarize  the  user  with  sta¬ 
tistical  results  and  to  present  some  possible  problems,  a  series  of 
graphs  comparing  observed  and  predicted  values  is  presented  in  Figure  3. 
The  first  graph  (A)  represents  one  case  of  perfect  prediction.  In  such 
cases  the  value  for  the  normalized  mean  error  is  0.0,  and  for  the  relia¬ 
bility  index  it  is  1.0;  for  the  paired  t-test  for  means,  the  computed  T 
is  undefined  since  the  denominator  in  the  equation  equals  zero. 

19.  Other  cases  are  possible  where  the  value  for  T  is  undefined 
but  where  predicted  and  observed  values  are  not  equal.  This  is  illus¬ 
trated  in  the  next  nine  graphs  (B-J)  in  Figure  3.  In  still  other  cases, 
as  shown  in  graphs  K  and  L,  a  computed  T  may  equal  0.0,  signifying 
that  the  means  are  the  same,  but  individual  predictions  are  not  the 
same  as  observed  values.  If  the  user  based  his  judgment  of  the  model 

solely  on  this  statistic,  he  would  not  draw  correct  conclusions. 

2 

20.  The  coefficients  a  ,  b  ,  and  r  for  the  linear  equation 
are  undefined  in  graphs  A-J  because  there  is  only  one  value  for  the 

x  axis  and  at  least  two  values  are  needed  for  the  computation.  For  the 

coefficients  for  the  linear  equation  to  represent  a  case  of  perfect  pre- 

2 

diction,  a  must  equal  0.0,  b  must  equal  1.0,  and  r  must  equal  1.0. 
It  is  possible  to  have  one  or  two  of  the  coefficients  correct  and  still 
have  predicted  values  which  are  not  equal  to  measured  values;  graph  l 
is  an  example  of  this  problem. 

21.  In  graphs  B-D  the  differences  between  the  predicted  and  ob¬ 
served  values  are  equal,  but  the  values  for  the  reliability  index  and 
the  normalized  mean  error  are  not  necessarily  the  same.  For  both  of 
these  statistics  the  values  are  scale-variant;  th8t  is,  the  same  numer¬ 
ical  difference  between  their  observed  and  predicted  values  will  produce 
a  smaller  calculated  value  as  the  numbers  being  coaq>ared  become  greater. 
In  the  case  of  the  normalized  mean  error,  scale  variance  depends  in  part 
on  the  observed  values,  since  they  occur  in  the  denominator;  thus,  the 
value  for  the  normalized  mean  error  is  the  same  in  graphs  B  and  D,  but 
different  from  the  value  in  C.  This  is  not  the  case  for  the  reliability 
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Value 

Figure  3.  Graphs  comparing  observed  (0)  and  predicted  (X)  values. 
Undefined  values  are  represented  by  U  (Sheet  1  of  4) 


Observed(O)  Predicted(X)  Reliability  Normalized 

Mean  Mean  T  Index  Mean  Error 


Linear  Equation 


Figure  3.  (Sheet  3  of  4) 


Observed (0)  Predicted(X)  Reliability  Normalized  — p--r  E9uatiop 


Figure  3.  (Sheet  4  of  4) 


index,  so  graphs  B  and  C  produce  the  same  result  which  is  different  froai 
D's.  Graphs  E-J  are  presented  to  show  the  results  of  different  statis¬ 
tics  when  observed  and  predicted  values  are  one,  two,  and  three  orders 
of  magnitude  apart.  Because  the  reliability  index  shows  an  increasing 
value  which  does  not  depend  on  whether  the  observed  or  predicted  value 
is  greater,  it  appears  to  be  the  best  statistic  for  aggregating  results. 

22.  But  to  use  these  tests  correctly,  a  number  of  assumptions 
Must  be  satisfied  (Sokal  and  Rohlf  1969).  First,  samples  should  be  col¬ 
lected  at  randoM.  This  is  rarely  done  in  collecting  water  quality  data, 
for  it  is  srare  important  to  be  able  to  gain  information,  for  purposes 
other  than  model  evaluation,  about  the  system  under  study  by  sampling 

at  set  intervals  of  time  and  at  uniform  distributions  through  space; 
in  addition,  the  majority  of  the  data  are  taken  during  daylight  hours. 
Second,  the  samples  should  be  independent  of  one  another.  This  assump¬ 
tion  is  violated  because  the  model's  prediction  at  one  point  in  time  de¬ 
pends  on  previous  values.  In  addition,  the  statistical  tests  are  valid 
only  if  the  samples  have  homogeneity  of  variances  and  are  normally  dis¬ 
tributed  and,  in  the  normal  case  of  regression,  if  the  independent  vari¬ 
able  is  measured  without  error. 

23.  Tests  and  possible  alternatives  for  these  assumptions  are 
available  (see  for  example,  Sokal  and  Rohlf  1969).  Some  of  the  adjust¬ 
ments  may  include  changing  the  basic  design  of  the  experimental  program, 
or  not  using  all  of  the  collected  data.  Rather  than  spending  time  and 
resources  in  checking  these  assumptions  or  deleting  data  to  make  sure 
that  the  statistical  results  are  absolutely  valid,  the  statistics  are 
here  used  as  a  tool  in  comparing  alternative  formulations.  A  listing  of 
the  statistical  package  and  information  concerning  its  use  are  included 
in  Appendix  B. 

24.  All  of  the  statistical  tests  which  are  proposed  compare  one 
predicted  value  against  one  observed  value,  but  most  reservoir  data  sets 
have  more  than  one  station  at  which  data  were  collected.  The  model  can 
be  run  stochastically,  a  process  which  produces  a  series  of  values  for 
each  variable  at  each  depth  for  each  time  step.  However,  the  cost  to 
run  the  model  stochastically  to  test  a  number  of  algorithms  for  the 


same  process,  or  to  add  or  delete  processes  or  variables,  would  be  pro¬ 
hibitive.  Also,  creating  a  value  that  is  supposed  to  be  measured, 
either  by  arithmetic  averaging  or  averaging  by  volume  weighting,  could 
produce  a  data  set  that  is  not  realistic.  For  example,  suppose  that 
the  water  entering  a  reservoir  is  colder  than  the  reservoir  proper. 

The  time  is  early  spring,  the  reservoir  is  well  mixed,  and  the  water 
is  entering  as  a  plug  flow,  warming  as  it  moves  through  the  reservoir. 
Because  the  reservoir  is  deeper  near  the  dam,  averaging  the  temperature 
at  different  stations  would  produce  a  data  set  that  has  cooler  surface 
water  over  a  warmer  hypolimnion.  Physically  this  does  not  happen,  and 
using  these  data  may  result  in  the  acceptance  of  inferior  algorithms. 

For  the  proposed  evaluation,  predicted  results  usually  come  from  run¬ 
ning  the  model  in  the  deterministic  mode,  and  comparisons  are  made  to 
data  collected  at  a  single  station,  usually  in  the  deepest  part  of  the 
reservoir.  Occasionally  the  model  will  be  run  in  a  stochastic  fashion, 
and  for  specific  variables  at  certain  depths  the  results  will  be 
graphed  with  measured  values  at  all  stations. 

Comparing  measured 

and  predicted  flux  values 

25.  Most  of  the  literature  concerning  the  comparison  of  model 
predictions  and  measured  values  deals  with  the  comparison  of  the  mass 
or  concentration  of  variables.  However,  it  is  possible  to  predict  the 
same  values  for  state  variables  with  very  different  fluxes  (the  term 
flux  refers  to  the  amount  of  materials  transferred  between  model  vari¬ 
ables;  e.g.,  the  amount  of  phytoplankton  ingested  by  zooplankton).  An 
example  of  this  is  shown  in  Figure  4.  In  each  model  there  are  four 
compartments,  with  the  values  for  the  initial  conditions  listed  above 
the  predicted  values  after  one  time  step.  For  both  models,  the  pre¬ 
dicted  values  were  the  same  even  though  the  fluxes  were  different.  TLe 
predicted  values  were  calculated  using  the  equations  listed  in  Table  2 
and  the  coefficients  and  driving  variables  listed  in  Table  3.  This  same 
type  of  problem  was  shown  by  Scavia  (1980)  to  occur  when  using  a  model 
of  Lake  Ontario.  Even  in  the  case  where  different  data  sets  are  used 
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Table  2 

Equations  Representing  the  Systems  Used  in  Figures  4  and  5 


dA 

dt 


’ACjVj  +  Dc4v4  -  Av6 


dB  A  D 

dt  -  AclVl  ■  Bc2V2  ’  Bc5V5  " 


dD 

dt 


jp 

-tt  ~  Bc0v0  -  Cc0v0  -  Cv, 
at  2  2  j  j  o 


Bc5V5  +  Cc3V3  '  Dc4v4  - 


+  I , 


Bv6  +  1 


+  1/ 


Dv6  +  1 


B 


D 


where 

t  equals  time 

A  ,  B  ,  C  ,  and  D  represent  mass  of  variables 
,  c4  ,  and  c,.  are  coefficients 
Vj  ,  ,  v4  ,  v,.  ,  and  v^  are  driving  variables 

1A  ,  lg  ,  1^,  ,  and  Ig  represent  the  mass  entering  respective 

compartments  from  outside  the  system 


Table  3 

Coefficients  and  Driving  Variables  Used  for  the 
Systems  Represented  in  Figures  4  and  5 


Coefficients 

Model  1 

Model  II 

Driving 

Variables 

Calibration 

Verification 

C1 

c2 

0.2429 

0.10 

V1 

1.0 

0.7368 

0.160 

0.2776 

V2 

0.7589 

0.506 

C3 

0.0157 

0.03114 

V3 

9.4 

10.12 

c4 

1.70 

0.5889 

v4 

0.1607 

0.5635 

c5 

5.70 

0.7440 

v5 

0.0360 

0.0239 

v6 

0.10 

0.091 

MODEL  I 


MODEL H 

Figure  5.  An  example  of  verifying  two  models 
having  the  same  initial  and  final  values  but 
different  flux  values 


for  calibration  and  verification,  there  is  no  way  of  telling  if  the  sec¬ 
ond  data  set  is  different  enough  to  allow  for  the  discovery  of  the  flux 
problem.  Suppose,  for  example,  that  a  second  data  set  was  available  to 
be  used  for  verification  (Figure  5).  Since  the  two  models  of  the  same 
system  had  different  coefficients  and  therefore  different  fluxes,  only 
one  of  the  models  would  be  expected  to  accurately  predict  the  measured 
values  for  the  verification  exercise.  The  verification  data  set  had 
different  driving  variables  (Table  3)  and  initial  conditions  (Figure  5) 


24 


from  the  calibration  data  set.  The  same  model,  outlined  in  Table  2,  and 
the  same  coefficients  (Table  3),  were  used  with  the  second  data  set.  As 
can  be  seen  in  Figure  5,  different  fluxes  were  predicted,  but  again  the 
two  models  predicted  the  same  values  after  a  single  time  step. 

26.  This  example  illustrates  that  one  cannot  depend  solely  on 
comparing  predicted  versus  measured  quantities  of  state  variables  when 
calibrating  or  verifying  models.  In  order  to  ensure  reliable  models  for 
a  particular  system,  calibration  and  verification  procedures  should  in¬ 
clude  comparisons  of  measured  versus  predicted  flux  values  as  well  as  of 
measured  versus  predicted  quantities  for  state  variables. 

27.  For  CE-QUAL-R1,  a  peripheral  package  has  been  developed  which 
provides  estimates  of  the  flux  values  used  to  calculate  values  for  vari¬ 
ables  for  the  next  time  step.  The  flux  values  reported  are  estimates; 
this  is  because  in  all  cases  for  the  flux  package  an  Euler  technique  is 
used  to  solve  the  equations,  whereas  most  variables  included  in  the  full 
water  quality  model  are  solved  using  a  2-step  (predictor-corrector) 

Euler  procedure.  The  connection  between  CE-QUAL-R1  and  the  flux  package 
is  a  file  created  during  the  simulation.  The  name  of  the  file  is  in¬ 
cluded  on  the  FILES  card  (Environmental  Laboratory  1982)  after  the  names 
for  the  plot  files.  Information  concerning  the  use  of  the  flux  package 
is  included  in  Appendix  C.  A  computer  listing  is  not  included  because 
the  computer  code  will  be  changed  as  new  algorithms  concerning  processes 
replace  present  formulations. 

28.  The  methods  for  model  evaluation  discussed  in  this  report  can 
be  used  when  the  model  CE-QUAL-R1  is  applied  to  either  an  existing  or  a 
proposed  reservoir.  In  the  case  of  a  preimpoundment  study,  however,  the 
data  for  use  in  model  calibration  and/or  verification  will  be  drawn  from 
other  nearby  reservoirs  judged  to  be  sufficiently  similar  to  the  pro¬ 
posed  project  that  observations  from  these  other  systems  are  represen¬ 
tative  of  the  reservoir  under  study.  Otherwise,  evaluation  of  model 
predictions  must  rely  strictly  on  scientific  and  engineering  judgment. 
The  reader  is  referred  to  the  reports  of  Fort  et  al.  (1977,  1979)  and 
Thornton  et  al.  (1976,  1977)  as  examples  of  uses  of  developmental 
versions  of  CE-QUAL-R1  in  p re impoundment  investigations. 
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APPENDIX  A:  STATISTICAL  EQUATIONS 


1.  Let  A.  equal  the  measured  value  at  the  n**1  depth  during 
th 

the  t  saaple  period  for  a  particular  variable,  and  let  P  equal 

the  value  for  the  corresponding  predicted  variable  where 

t  =  1,...,T  (total  nuaber  of  sampling  periods  for  a  variable) 

n  =  1,...,N  (total  nuaber  of  depths  saapled  for  each  saapling 
period) 

N  =  total  nuaber  of  observed  samples  for  all  depths  for  all 
saapling  periods 

2.  The  reliability  index  RI  (Leggett  and  Williams  1981)  is 
defined  as: 

1 

RIt  =  - 
1 


for  a  sampling  period  or 

1  + 

RI  =  - 

1  - 

for  all  sampling  periods  for  a  particular  variable. 


A1 


The  paired  t  test  for  means  B  is  calculated  as 


for  a  sampling  period  or 


for  all  sampling  periods  for  a  particular  variable. 

3.  The  normalized  mean  error  NME  (Gordon  1981;  see  also 
Wlosinski  1982)  is  calculated  as 


for  all  sampling  periods  for  a  particular  variable. 


A2 


4.  The  final  statistical  test  is  a  regression  analysis  of  the 
relationship  between  predicted  and  measured  values.  The  statistical 
test  gives  estimates  of  a  and  b  in  the  equation 


a  +  bA. 

tn 


which  could  then  be  compared  with  the  equation  for  perfect  prediction, 

where  a  equals  0.0  and  b  equals  1.0,  using  the  Students  t  distri- 

^  2 
bution.  In  addition,  the  square  of  the  correlation  coefficient  (r  ), 

which  is  a  measure  of  the  variance  accounted  for  between  observed  and 

predicted  values,  can  be  calculated. 

5.  Calculate 


T 
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The  coefficients  b 


Further  calculate 


T  -  M 


A  = 


Si 


p  =  i? 

r  IN 


IA 


FA  = 


(1AY 

~sr 


(INT  -  1 


5p2  .  OP) 

FP  -  » 

FP  ■  Tar- 1" 


ZAP  - 


FAP  = 
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IAIP 

"sr 


and  a  are  calculated  from  the  equations 


b 


FAP 

FA 


a  =  P  -  bA 


SD  =  [(IN)  -  1 ]b2FA 

ST  =  [(IN)  -  1]FP 
„  _  ST  -  SD 

S  -  TINT-- 2 
^  =  [ (IN)S-  1 ]FA 


A4 


Va 


SHUN)  -  1]FA  +  A2IN) 
IN[(IN)-  1 ]FA 


The  null  hypothesis  NH  ,  that  the  computed  slope  equals  one,  is  tested 
with 


NH  = 


b  -  1 
■^Vb 


with  (IN)  -  2  degrees  of  freedom.  The  null  hypothesis,  that  the 
intercept  equals  zero,  is  tested  with 

NH  = 

\Va 


2 

with  (IN)  -  2  degrees  of  freedom.  The  correlation  coefficient  r 
equals 
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APPENDIX  B:  STATISTICAL  PACKAGE 


1.  As  with  CE-QUAL-R1,  the  statistical  package  resides  on  the 
Boeing  Mainstream-EKS  interactive  time-sharing  computer  system.  A  user 
must  have  his  own  account  number  and  should  have  some  knowledge  of 
Boeing  procedures.  Because  of  the  large  amount  of  data  generated  by 
the  CE-QUAL-R1  model,  two  executions  are  needed.  The  first  execution 
is  used  only  to  delete  predicted  data  on  dates  when  no  measured  data 
are  available.  To  execute  this  program,  prepare  a  file  as  follows  for 
use  by  the  SUBMIT  directive: 

GRAXEQ , CM320000 , T300 , PI . 

USER , ID , PASWORD . 

GET , GRA0BJ/UN=CER0B5 . 

GET , TAPE5=LDATE . 

GET , TAPE89=PLDG14 . 

FILE ,TAPE6 ,FF=YES . 

LOADXEQ , F=GRAOBJ , M=FULL/MAPGRA . 

EXIT.U. 

REPLACE, MAPGRA. 

EXIT,U. 

REPLACE  ,TAPF.6=GRADAT . 

EXIT.U. 

REPLACE , 0UTPUT=PUT0UT1 . 

EXIT.U. 

COST , LO=F . 

EXIT.U. 

DAYFILE.DFXGRA. 

REPLACE, DFXGRA. 

2.  Names  that  are  underlined  can  be  changed  at  the  users  discre¬ 
tion.  File  LDATE  contains  the  dates  on  which  measured  data  are  avail¬ 
able.  One  date  should  be  in  the  first  six  columns  of  each  line:  the 
first  two  columns  represent  the  last  two  digits  of  the  year,  columns 
three  and  four  contain  the  number  of  the  month,  and  columns  five  and 


six  contain  the  day  of  the  month.  File  PLDG14  is  the  same  file  as  is 
used  by  the  interactive  graphics  package  (see  Environmental  Labora¬ 
tory  1982). 

3.  Upon  satisfactory  completion  of  the  execution,  four  files 
will  be  created  and  permanently  stored.  File  MAPGRA  contains  the 
storage  location  map.  File  GRADAT  is  the  output  file  that  will  be  used 
in  the  next  execution  which  will  perform  the  statistical  analysis.  File 
PUT0UT1  contains  information  concerning  problems  during  execution.  File 
DFXGRA  is  the  dayfile  which  describes  the  execution. 

4.  To  execute  the  statistical  package,  the  following  file  is 
needed. 

STSTXEQ,CM320000 ,T15 ,P01 . 

USER , ID , PAS WORD . 

GET , STSTOBJ /UN=CER0B5 . 

GET,TAPE22=VD794. 

GET , TAPE23=GRADAT . 

GET , TAPE6=STSWICH . 

FILE ,TAPE7 ,FF=YES . 

LOADXEQ , F=STSTOBJ , M=FULL/ MAPSTST . 

EXIT,U. 

REPLACE, MAPSTST. 

EXIT,U. 

REPLACE ,TAPE7=STSTDAT . 

EX1T.U. 

REPLACE , OUTPUT=PUTOUT 1 . 

EXIT.U. 

COST,LO=F. 

EX1T,U. 

DAYF I LE , DFXSTST . 

REPLACE, DFXSTST. 

5.  File  VD794  is  the  file  that  contains  measured  values.  An  ex¬ 
ample  of  this  file  is  given  in  Figure  B1 .  The  first  40  characters  of 
the  first  line  are  for  informational  purposes;  the  rest  of  the  file  con¬ 
tains  measured  data.  All  of  the  data  for  a  particular  variable  should 
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Figure  Bl.  An  example  of  the  data  set  File  VD794  which  contains 
measured  values  used  by  the  statistical  package 


be  grouped  together;  within  each  group,  data  should  be  ordered  according 
to  date.  For  each  date  the  following  information  is  needed.  The  first 
line  contains  three  variables  describing  a  block  of  data.  Columns  5  and 
6  contain  a  code  number  for  each  variable.  (A  list  of  the  variables  and 
their  code  numbers  are  given  in  Table  Bl . )  Columns  8  through  13  contain 
the  date  on  which  the  block  of  data  was  measured:  columns  8  and  9  con¬ 
tain  the  last  digits  of  the  year,  columns  10  and  11  contain  the  numeri¬ 
cal  description  of  the  month,  and  columns  12  and  13  contain  the  day  of 
the  month.  Columns  15  through  17  contain  the  number  of  data  points  for 
that  block  of  data.  Each  data  point  consists  of  a  pair  of  numbers:  the 
first  is  the  depth,  in  meters,  where  the  sample  was  measured;  the  second 
number  is  the  concentration  of  the  variable.  The  depth  is  measured  from 
the  surface,  and  the  data  are  ordered  from  the  bottom  to  the  surface. 
There  are  four  pairs  of  numbers  on  each  line  in  columns  19-22  and  24-32, 
34-37  and  39-47,  49-52  and  54-62,  and  64-67  and  69-77. 

6.  File  GRADAT  is  the  file  that  was  created  in  the  previous  exe¬ 
cution.  File  STSWICH  is  a  one-line  file  containing  information  on  vari¬ 
ables  for  which  statistics  will  be  performed.  A  in  a  particular 
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Table  B1 

Variables  and  Their  Codes  for  the 
Statistical  Package 


Code 

Variable 

1 

Temperature 

2 

Zooplankton 

3 

Algae  1 

4 

Algae  2 

6 

Total  manganese 

7 

Detri tus 

8 

Dissolved  organic 

matter 

9 

Ortho-Phosphate  - 

P 

10 

Inorganic  carbon 

11 

Ammonia  -  N 

12 

Nitrite  -  N 

13 

Nitrate  -  N 

14 

Oxygen 

15 

Carbon  dioxide 

16 

PH 

17 

Alkalinity 

18 

Total  dissolved  solids 

19 

Suspended  solids 

20 

Total  iron 

21 

Sulfate 

22 

Reduced  manganese 

23 

Reduced  iron 

24 

Iron  sulfide 

25 

Reduced  sulfide 

26 

Coli forms 

column  will  cause  the  statistical  analysis  to  be  performed.  The  column 
code  for  variables  is  presented  in  Table  B1 . 

7.  Upon  successful  completion  of  the  execution,  four  files  will 
be  created  and  permanently  stored.  File  MAPSTST  contains  the  storage 
location  map.  File  STSTDAT  contains  statistical  output  of  which  an  ex¬ 
ample  is  given  in  Figure  B2.  File  PUT0UT1  contains  information  con¬ 
cerning  problems  during  execution.  File  DFXSTST  is  the  dayfile  which 
describes  the  execution. 

8.  A  listing  of  the  statistical  package  is  presented  in 
Figure  B3. 


BA 


STATISTICAL  PR06RM  FOR  CE-0UW.-R1  COMPARING  FILE 


Figure  B2.  An  example  of  output  from  the  statistical  package 


PROGRAM  STSTIC(TAPE22. TAPE23.TAPE6. TAPE7) 

STATISTICAL  PROGRAM  FOR  CE-QUAL-R1  COMPARING  MEASURED 
AND  PREDICTED  VALUES 

DIMENSION  niTLE(2»4).K0DE0(450).KDATQ(45O>»KPT0(450). 

*  DEP0(450.50>»VALU0(450.50).NAMES<26) »NSWICH(50> ,VALU(70) » 

*  VALUP(70) .DEPP (70) 

DATA  NAMES/* TEMP* » 1  ZOOPLANK 1 » ’ALGAE! * » * ALGAE2* , * ALGAE3* . 

*  *  MN2+MN4* .’DETRITUS* t'D.O.M.’t*  P04-P* » ’ T . I «C .  * » ’NH4-N* » 

*  ’ N02-N’. ’ N03-N ’. ’OXYGEN*. ’C02*.’ PH’. ’ALKLINITY* . ’T.D.S. ’ . 

*  ’ S. S.l’.’T, IRON’.  ’SULFATE ’ . *R.MN’ . ’R.FE’ . ’FE-S’ . ’SULFIDE’. 

*  ’COLIFORM’/ 

DATA  KODOLD/l/i NS TOP/O/ 

READ (22. 30) (ITITLE(1»I).I=1.4) 

Rw®>?»imMKiusm!n).(iTITLE(2.I), 1-1.4, 

29  FORMAT ( 1H1 » 40X. 

450HSTATISTICAL  PROGRAM  FOR  CE-QUAL-R1  COMPARING  FILF  >•. 

*  / . 20X . 4A 1 0 . 3X . 3HAND , 3X . 4 A 1 0 . / / / ) 

HRITE( 7.27) 

27  FORMAT (9X.39HRELIABILITY  INDEX  LEGGETT  AND  MILL  I AMS 

*  . 37H  ECOLOGICAL  MODELLING  13(1981)303-312./. 

49X.53HPAIRED  T  TEST  FOR  MEANS  SOKAL  AND  ROHLF  1969  BIOMETRY. 
*24H  kl.H. FREEMAN  AND  COMPANY./. 

49X.21HN0RMALIZED  MEAN  ERROR. 

*25H  (SUM(AB3(P-0)/0)*100. ,/N. 

*/»  9X»  47HC0EFFICIENTS  FOR  THE  LINEAR  REGRESSION  EQUATION.///) 
READ (6. 28) (NSWICH ( J) . J=1 ,50) 

28  FORMAT (5011 ) 

30  F0RMAT(4A10) 

1  =  0 

32  1=1+1 

REAIK  22. 35)K0DE0( I ) *  KDATO ( I )»KPTO(I).KT» 

*  <DEP0(I.K),VALU0(I,K).K=1.KT> 

IF (EOF ( 22 ) . NE . 0)G0  TO  34 

IF (KT.EQ.4 ) BACKSPACE  22 
GO  TO  32 

34  11=1-1 

35  FORMAT (4X.I2.1X.I6.1X.I3.T15.I3.1X.4(F4»1.1X.F9>3>1X)./. 

*  (18X.F4.1»lX.F9.3.1X»F4.1.1X.F9.3.1X.F4.1,lX»F9.3.1Xi 

*  F4.1. 1X.F9.3) ) 

40  READ(23»35)K0DEP»KDATP.KPTP»KT » ( DEPP < I ) , VALUP( I )  *  1  =  1 »KT  > 
IF(E0F(23) *NE.O)GO  TO  198 
IF (KT .EQ. 4 ) BACKSPACE  23 
IF(KODOLD.NE.KODEP)GO  TO  200 
44  IF(NSUICH(KODEP) . NE,0)G0  TO  46 
KODOLD=KQDEP 
GO  TO  40 
46  DO  55  1=1.11 
JJ=I 

IF(KODEP.NE .KODEO( I ) .OR. KDATP.NE. KDATO ( I ) > GO  TO  55 
GO  TO  74 
55  CONTINUE 
GO  TO  40 
74  J=0 

FIGURE  WHICH  LAYER  PREDICTIONS  WILL  BE  USED 


KT=KPTO( JJ) 


Figure  B3. 


A  listing  of  the  statistical  package 
(Sheet  1  of  4) 
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DO  80  I=i f KT 
76  J-J+l 

IF(DEPP< J)  .  GT . DEPO<JJ>  I )  )G0  TO  76 
DI=ABS\DEPP( J ) -DEPO( JJ* I ) ) 

IF<ABS<BEPP< J-1)-DEP0(JJ.I)).LT.DI)J=J-1 
VALU< I)=VALUP< J) 

J=J-1 
80  CONTINUE 

* 

•  T  STATISTIC  FOR  MEANS 

* 

RDAL2=0. 

RBAL=0  4 
0B=  0 
PR-0 

HO  82  I=1»KT 

RDAL2=RDAL2t(VALU0(JJ,I)-VALU<I>>**2. 

RDAL=RDAL+VALUO( JJtl)-VALU(I) 

OB=OB+VALUO< JJ.I) 

82  PR=PR+VALU< I ) 

OBT=OBT+OB 
PRT=PRT +PR 
OBA=OB/KT 
r'RA=PR/KT 
KTTOT=KTTOT+KT 
RS2T=RS2T+RDAL2 
RST=RST+RDAL 
ZT=KT 

ZRT=<0BA-PRA)/SGRT(<(ZT*RDAL2>-RDAL**2.)/<ZT**2.*<ZT-1.))> 

S=0 

SK=0 

* 

*  RELIABILITY  INDEX 

* 

DO  90  1-1 »KT 

IF(VALU<I).EQ.0)G0  TO  90 
SK=SN+1. 

R=VALUO( JJi 1 5 /VALU< I ) 

S=((l.-R)/(l.+R))**2fS 
90  CONTINUE 
STOT=STOT+S 
SKTOT=SKTOT+SK 

IF(SK.EG.O)GO  TO  96 
AB=1.+<SQRT<<1,/SK)*S>> 

AC=1 . -(SORT <  < 1 . /SK )  *S) ) 

AK=AB/AC 

* 

*  NORMALIZED  MEAN  ERROR 

* 

96  P=0 
SL=0 

DO  120  1=1 rKT 

I F ( VALUO ( J J » I > .EQ. 0 )G0  TO  110 
SL*SL+1 . 

P=(ABS(VALU<I)-VALUO(JJrI))/VALUO(JJ»I))+P 
110  CONTINUE 
PTOT=PTOT+P 
SLTGT=SLTOT+SL 

IF(SL.EQ.O)GO  TO  111 
F=(P*100. )/SL 
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******** 


* 

* 

* 

* 


COEFFICIENTS  FOR  THE  LINEAR  EOUATION 
CALCULATED  ONLY  OVER  ALL  TINES 

DO  300  1=1 »KT 
SUMA=SUMA4VALU0( JJ, I ) 

SUMA2=SUMA2+VALU0( JJ, I )**2. 

SUMP=SUMP+VALU<I) 

SUMP2=SUMP2+VALU( I ) **2 , 

SUMAP=SUM AP+VALUO « J J» I ) *VALU  < I ) 

300  CONTINUE 
SUMKT=SUMNT+KT 
IF(NSNCK2.£Q.0)URITE(7f 114) 

NSWCH2=1 

114  FORMAT ( 1H0 > 2X 1 8HVARI ABLE > 5X > 4HDATE » 5X » 9HNUNBER  OF . 4X»  8H0BSERVED. 

*  5X,9HPREDICTED,7X,1HT,7X,11HRELIABILITY,3X,10HN0RMALIZED,/,14X. 

*  6HYYMMDD,3X»llHC0MPARIS0NS»5X»4HMEAN»10X»4HMEAN»6Xf 9HSTATISTIC 

*  f 5X»5H INDEX »7X» 10HMEAN  ERROR) 

111  URITE<7 , 112 )NAMES(KODEP) .KDATP .KT.OBA.PRA.ZRT.AK.F 

112  FORMAT (lXfA10t3X»I6f6XfI4»AXr5<E10.3f3X)) 

CALCULATE  STATISTICS  FOR  A  PARTICULAR  VARIABLE 
OVER  ALL  TIMES  AND  DEPTHS 

GO  TO  40 
198  NST0P=1 
200  CONTINUE 

IF<NSUICH(KODOLD) . NE.DGO  TO  214 

A=0BT/KTTO1 

B=PRT/KTTOT 

AB= 1 .  +  ( SORT ( ( 1 . /SKTOT ) *STOT ) ) 

AC=1  .-(SORT <(l./SKTOT)*STO'i  >> 

AK=AB/AC 

F=<PT0T*100.)/SLT0T 

ZRT=<A-B)/SORT< ( <SKT0T*RS2T)-RST**2, )/ 

*  (SKT0TM2.XSKT0T-1.))) 

AVEA=SUMA/SUNKT 

AVEP=SUMP/SUMKT 

FA= ( SUMA2-SUMAM2 , /SUMKT ) / ( SUNKT-1 .  > 
FP=(SUMP2-SUMP**2./SUMKT)/(SUMKT-  1 . ) 

FAP= ( SUMAP-SUMAtSUMP/SUMKT  >/< SUMKT- 1 . ) 

ESTB=FAP/FA 
ESTA=AV£P-ESTB*AVEA 
SD=<SUMKT-l.)*ESTB*t2.*FA 
STXSUMKT-l.  >*FP 
SI=(ST-SD)/<SUMKT-2« ) 

VB=SI/< ( SUMKT- 1 . )*FA) 

VA=SI*( (SUMKT-1 . )*FA+SUMKTIAVEA*»2, )/ 

*  < SUMKT X SUMKT- 1 . )*FA) 

TSL0PE=<ESTB-1 * )/SQRT  <  MB ) 

TNTCPT=ESTA/SQRT (VA) 

RSQ=SD/ST 

STDRD2=  ( SUMP2-SUMPM2 ,  /SUMKT )  -ESTBM2 .  *  <SUNA2-SUMA**2 .  /SUMKT ) 
STDRD=S0RT(STDRD2) 

URITE(7. 121 ) NAMES (KODOLD) .KTTOT  . A.B.ZRT ,  AK.F 
121  FORMAT ( 1H0. A10. 10H  ALL  DATES, 5X» I4,6Xr5<E10.3, 3X) ) 

WRITE<7,301 )ESTA,ESTB,TSLOPE,TNTCPT  .RSQiSTDRD 

301  FORMAT ( IX , 2HA= » E9  *  3 , 2X » 2HB= , E9 . 3 , 1 5H  T  FOR  SLOPES®, E9.3, 

*  20H  T  FOR  INTERCEPTS® »E9, 3,1 1H  R  SQUARE®, E9, 3, 

*  19H  ST, ERROR  OF  EST.®,E9.3> 
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1F(NSTQP,EQ.1)G0  TO  500 

NSWCH2=0 

0BT=0. 

PRT=0. 

KTTOT'O. 

SKTOT-O. 

STOT:rO. 

SLT0T=0. 

RS2T"0» 

RST=0. 

PT0T=0. 

SUNA=0. 

SUMA2=0. 

SUMP=0. 

SUMP2=0. 

SUHAP=0 • 

SIMKT^O, 

214  KODOLD---KGBEP 
60  TO  H 
500  STOP 
END 
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APPENDIX  C:  FLUX  PACKAGE 


1.  As  with  CE-QUAL-R1,  the  flux  package  resides  on  the  Boeing 
Mainstream-EKS  interactive  time-sharing  computer  system.  A  user  must 
have  his  own  account  number  and  should  have  some  knowledge  of  Boeing 
procedures.  To  execute  the  flux  package,  prepare  a  file  as  follows 
for  use  by  the  SUBMIT  directive: 

FLXXEQ , CM300000 ,T300,P01. 

USER, ID, PAS WORD. 

GET , FLXOB J/UN=CER0B5 . 

GET,TAPE61=PLDGFLX. 

GET , TAPE5=FLXDDAY . 

FILE, TAPE6 , FF=YES . 

LOADXEQ , F=FLXOBJ ,M=FULL/MAP2 . 

EXIT,U. 

REPLACE, MAP2. 

REPLACE , TAPE6=0UTPF . 

REPLACE , OUTPUT=PUTOUT 1 . 

EXIT.U. 

COST , LO=F . 

EXIT.U. 

DAYF I LE , DAYFX . 

REPLACE, DAYFX. 

2.  Names  that  are  underlined  can  be  changed  at  the  user's  discre¬ 
tion.  File  PLDGFLX  is  the  file  on  which  information  from  a  CE-QUAL-R1 
simulation  is  stored.  File  FLXDDAY  (see  Figure  Cl)  is  a  four-line  file 
containing  information  needed  for  output.  The  first  line  contains  the 
hour  of  the  year  for  which  output  is  requested;  up  to  16  values  can  be 
specified  in  fields  of  five  characters  each.  The  second  line  contains 
either  a  blank  or  a  '1'  in  the  first  16  columns.  A  '1'  signifies  that 
information  is  needed  for  a  particular  variable;  the  16  variables  for 
which  information  is  gathered  are  listed  in  Table  Cl. 


Cl 


1  1  2 
Column  15050 

1032  1056  1296  1536 
111  111111 
SURFACE  CUBICM 
1  168  24 

Figure  Cl.  An  example  of 
file  FLXDDAY 


Table  Cl 


Variables 

for  Which  Flux  Information 

is  Available 

Column 

Variable 

1 

Fish  1 

2 

Fish  2 

3 

Fish  3 

4 

Benthos 

5 

Zooplankton 

6 

Algae  1 

7 

Algae  2 

8 

Detritus 

9 

Sediment 

10 

Dissolved  Organic  Matter 

11 

Ortho-Phosphate  -  P 

12 

Carbon 

13 

Ammonia  -  N 

14 

Nitrite  -  N 

15 

Nitrate  -  N 

16 

Oxygen 

3.  The  third  line  of  FLXDDAY  contains  two  variables.  The  first 
variable  concerns  how  the  fluxes  are  to  be  summed  according  to  layers. 
(It  must  be  remembered  that  CE-QUAL-R1  contains  a  variable-layer  scheme 
and  that  the  layers  are  numbered  from  the  bottom  to  the  surface. 
Throughout  the  year  the  number  of  layers  may  change,  so  that  the  sur¬ 
face  layer  may  not  always  have  the  same  layer  number.  If  the  user  is 
interested  in  a  particular  process  occurring  at  the  surface — for  exam¬ 
ple,  surface  aeration — he  would  have  to  look  at  different  layer  numbers 


during  different  times  of  the  year.  To  make  it  easier  to  study  the 
fluxes  occurring  in  the  epilimnion,  it  is  possible  to  sum  fluxes  in  re¬ 
lation  to  the  surface  by  putting  the  word  SURFACE  in  columns  1  through  7 
of  the  third  card.  Otherwise,  the  fluxes  will  be  summed  from  the  bottom 
of  the  reservoir).  The  second  variable  concerns  units  of  output:  if 
the  word  CUBICM  is  specified  in  columns  11  through  16,  fluxes  will  be 
reported  in  g/cu  m;  fluxes  will  be  reported  in  units  of  kg/layer  with 
any  other  specification. 

4.  The  fourth  line  of  FLXDDAY  contains  three  variables,  the  first 
of  which  concerns  the  number  of  time  steps  for  which  fluxes  are  to  be 
accumulated;  the  value  is  shown  in  columns  1  through  5,  right- justified. 
For  example,  if  the  original  simulation  used  a  24-hour  time  step,  a  '1' 
in  column  5  would  give  information  on  fluxes  on  a  daily  basis;  a  value 
of  30  would  accumulate  fluxes  for  periods  of  approximately  1  month.  The 
second  variable,  in  columns  6  through  10,  represents  hours  and  allows 
output  to  be  specified  at  particular  intervals  (this  information  is  out¬ 
put  in  addition  to  that  specified  in  the  first  line).  The  third  vari¬ 
able,  in  columns  11  through  15,  specifies  the  Julian  date  for  the  start 
of  the  simulation. 

5.  Upon  satisfactory  completion  of  the  execution,  four  files  will 
be  created  and  permanently  stored.  File  MAP2  contains  the  storage  loca¬ 
tion  map.  File  OUTPF  is  the  normal  output  from  the  flux  package;  an 
example  of  output  is  given  in  Figure  C2.  File  PUT0UT1  contains  informa¬ 
tion  concerning  problems  during  execution.  File  OAYFX  is  the  dayfile 
which  describes  the  execution. 


C3 


SIMULATION  HOUR-  744  LOWEST  NUMBER  OF  LAYERS1  55  LAYERS  NUMBERED  FROM  THE  SURFACE 
FUKES  ARE  IN  KG/LAYER  ACCUMULATED  OVER  THE  LAST  7  SIMULATION  TIME  STEPS 

mmnttmmmmmmuummnmtmmtmmmmmmmntmm 

ttmummmmmntmmmmmtmmmttmmmmtnmmmutu 


“8>  * 


|  I  I  I  I  I  I  I  I  I  I  I  I  t  I  I  I  I  4  I  I  I  I  I  I  I  I  I  I  I  I  I 


MrtrOrtMMMM  MWMMMWPOMiMMMWMPOWrtMrOMMMrOMMM 

SOOOOOOOO OOOOOOOOOOOOOOOOOOOOOOOOO 
+  4  ■4"4»4'  44  4  4  4  4  4  4  +  *►  4-  "t”  -*>  +  ■+■  4  ^  *♦^7'. 


a  85 


t  i  i  i  i  i  i  i  i  ■  i  i  i  i  >  i  i  <  ■  i  i  i  i  i  i  i  t  i  i  i  i  i  • 


m 

<4- 

c 


CW  fO  <M  CW  CW  CW  CWCW  CW  CM  CWCW  CM  Csl  CWCN  CM  CM  CWCWCWCW  CM  CM  CWCW  CWCWCWCW  CWCWP* 
zzoo oooooo oooooooooooooooooooopoooo 
CD  2  ♦  4*  4*  *4*  4  4*  ^  41 4-  -4-  •4“  -4“  4-  -4-  4  4*  41  4  4  *4-  4“  4  4-  -4-  4-  4  4-  4 4-  "4- 

3  8>  I*  I  t  I  I  I  »  I  •  »  I  I  •  I  *  «  I  I  I  4  I  I  I  I  •  t  I  I  I  4  I  I  * 

OU4 


a 

5 


X4444444 

OOOOOOOO 
►4  -4  4-  4-  4-  -4  4-  4 


ooooooooooo< 


>000000000 


81 


I  I  I  I  I  I  4  I  I  I  I  I  <  t  I  I  •  I  I  I  I  I  I  I  I  I  I  I  I  |  !  I  I  I 


a 


S4  •— I  •*— «  yS  4  4 44 44 ^«44 4444 44444 44444 4 4  O  OOP 

OOOOOOOO oooooooooooooooooooooooop 

4  4-  4“ 4- 4- 4- 4- 4- 4-  4- 4 4 4 4 4 4 4 4 4 4 4 4 4  4 4 4 4 4 4 4 4 4 4 4 

£3K£!2r!8£p;&S25Ks1skP;f:BK§5K5SK222sS!S 

-♦miPSLrSiTW^-^  4F5ror3»rprocwrMcwcMcwcww^.- 


a  as 


«  i  i  i  >  >  i  >  »  i  i  »  »  i  *  J  i  i  r  >  i  j  i  i  i  i  i  i  i  i  r  i  i 


3 


CM  CW  CWCW  CWCW  CWCW 

OOOOOOOO  o o oo O O OO OOOOOO OOOOOOOO ooO 
44  4444  44  4 4444 44444444444444444 444 

•— irors  cm 4  'Cooorv  incocNangoo5ro»-irvroorxL‘>'OOOfNCO>owToo*HuTor^ 
>-oeiofv  rv'CinN-  4ro  cowcoroop<o*tn»Hrs»Ofv^^P‘'ON-^ 

^  o-  o- coajrvf^i «ouTinmT ■r^nrOMnncvtN 

co  ►-«  *  -  -  . 

3*  . 


I  I  I  I  <  t  I  I  I  I  (  I  I  I  I  t  I 


t  I  I  I  I  I  <  I  t  t  I  I 


i  .  <jr  Ul 1  **  UJ UJ  i  «,»■«!  i»ii ■  I  UJLWUJIU  Lii  UJ  \  V LaJI  Uj 1  ■ » UlliJUlUIUiLlJUJujuiUUlJUi  UJUJUJ 

a in c- 0  4 coo  '5in'OCp*-iM3^'r4*r  ^afs.^MSrnSrt^'orx'O^in 

'fl(NcsrifvcNp»o^  roYNmuirs 

IrtPrtOPMO'-O  MOfMrt(N©NlAWrAO‘NlF)MIMOfll)rv>O^WM^OM 

>ors  rv(S'0>c  u'Jin  w)UT^«f*r  •»  ro  ro  ro  r*o  cm  cm  cm  cm  cm  r  w 

.  . . . . . . . 

_  i  i  t  t  i  t  i  i  i  i  i  i  i  t  i  <  i  i  «  i  i  <  t  t  i  i  t  i  i  i  i  i  i 

44444444  44444  MMKlMWrtr»IWMtOWMMrOWWM  WMM 
ZOO OOOOOO OOOOOOOOOOOOOOOOOOOOOOOOO 
—  *■  -*-  *“444444444444444444444444444444 

. 

0-00 

—  o 

S3 


to  I  I 

3 


i  I  t  I  I 


I  I  I  I  I  I  I  I  l  I  I  l  •  l  I  l  I  • 


2SSSSSSSSSS2 


'8 


liOCO  r-v  -Cl  •<■!*>  CM  r-t>OW 


OOO oooooo o OOO OOO OOO OoO 


85 


-I  CM  m  4  UD  -o  r^oo  ©“  o 


o«HNP<xin'Orsoopor4C4to^-ii')PNmo>o«HMrii 
.,•*  ^^w^*-**-*^*-*^-' cwcwrwcw  rwcwrwrwcwrwrororor** 


C4 
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